/*
Master File for 'Undue Burden Beyond Texas'
Authors: Joanna Venator and Jason Fletcher
Date: 5/1/2021

This file replicates the analyses included in the paper, excluding analyses of births
which require restricted access NVSS data. Code to replicate those files is included, but commented out. For researchers who get access to the NVSS files and would like the code to create the .dta that we use in our analysis, email venator@wisc.edu to receive the do-files to clean the NVSS files.


*/

*Note: Comment out and comment in where ever you want files to end up saved when you forget to specify a specific path
// Add in your own directories (do not delete existing ones) and comment out directories you're not using with **
/*

*/


global home "[]\ReplicationFiles_VenatorFletcher_JPAM"


global logsave "${home}/logfiles"
global dtasave "${home}/dtafiles"
global rawdata "${home}/rawdata"
global WIdata "${home}/Wisconsin"

local restrict = 0 // set equal to 1 if using data with births; set equal to 0 if only using public use data

cd "${home}"

// Stata commands that may need to be installed if you don't already have them:
ssc install shp2dta, replace // package for shape files for map
ssc install spmap, replace // mapping package
ssc install estout, replace //  package for outputting files into csv/ tex files
 ssc install synth


// FIGURE 1

shp2dta using "${rawdata}/tl_2010_55_county10.shp", data("${rawdata}/WI_data")  coor("${rawdata}/WI_coordinates") genid(id) replace
use "${dtasave}/ClinicDataSet.dta", clear
egen dist2009 = mean(closest_dist_)if monthyear == 1, by(county) 
egen dist2017  = mean(closest_dist_)if monthyear == 102, by(county)
cap drop _merge 
cap drop changedist
preserve 
collapse (min) dist2017 dist2009, by (county)

gen changedist = dist2017- dist2009
tempfile temp1
 save "`temp1'" 
 restore
 merge m:1 county using "`temp1'" 
 preserve
 use "${rawdata}/WI_data", clear
 destring, replace
  g county = STATEFP10*1000+ COUNTYFP10
  
  tempfile temp2
   save "`temp2'" 
 restore
  

   preserve 
  import excel using "${rawdata}/WI_clinic_pointdata.xlsx", clear
  list


  rename A xcoord
  rename B ycoord
  rename C type
  label define type_l 1 "Open 2009-present" 2 "Closed During Analysis Period" 3 "Opened Post-Analysis Period" 4  "Opened During Analysis Period"
  label values type type_l
    save  temp3.dta, replace
  restore
 drop _merge
  merge m:1 county using "`temp2'" 
  preserve
  keep if monthyear == 1
 spmap changedist using "${rawdata}/WI_coordinates", id(id) fcolor(Blues) clmethod(custom) clbreaks(-56.7 -0.001 0.001 25 50 100 130.4)  ///
 legtitle("Change in Distance (miles)") title("Change in Distance to Nearest Abortion Clinic, 2009 to 2017", size(*0.8)) scalebar(units(1) scale(1) label(100 miles)) ///
    point(data("temp3.dta") xcoord(ycoord)               ///
    ycoord(xcoord)    by(type)        ///
    psize(absolute) fcolor(green red yellow black) shape(o X S T) size(*1.5 *2 *.5 *.75) legenda(on) legtitle("Clinic Type") leglabel("Open 2009-present" "Closed During Analysis Period" "Opened Post-Analysis Period" "Opened During Analysis Period") ) 
 graph export "WI_changedistancemap.png", replace
restore
//FIGURE 2
use "${dtasave}/ClinicDataSet.dta", clear
xtset county monthyear

local yearloc = 2009
local monthloc = 1
foreach num of numlist 1/134{
foreach n of numlist 13 25 37 49 61 73 85 97 109 {
if `n' == `num'{
local yearloc = `yearloc '+ 1 
}
}
local monthloc = `num' - 12*(`yearloc' - 2009)
label define monthyearlabel `num' "`monthloc'/`yearloc'", add

}
*label values monthyear_birth monthyearlabel
egen meandist = sum(closest_dist_), by(monthyear)
egen tag = tag(monthyear)
drop if monthyear >108

svyset county  [pw=pop]

svy: mean closest_dist_, over(monthyear)
sort monthyear
by monthyear: gen ET = sum(closest_dist_*pop) / sum(pop)
by monthyear: replace ET=ET[_N]
drop if monthyear == 0
//Average distance to nearest abortion clinic
twoway line ET monthyear if tag, sort ///
xlabel(1 "1/2009" 13 "1/2010" 25 "1/2011" 37 "1/2012" 49 "1/2013" 61 "1/2014" 73 "1/2015" 85 "1/2016" 97 "1/2017"  ) ///
scheme(s1mono) name("graph_meandist", replace) xtitle("Month-Year") ytitle("Mean Distance to Closest Abortion Clinic (mi)") ///
xline(56 82, lcolor(red))title("Residents' Average Distance to Closest Abortion Clinic")
graph export "Figure2_PanelA.png", replace


svyset county  [pw=pop]
g servicepop= servicepop_measure1/1000
svy: mean servicepop, over(monthyear)
sort monthyear
by monthyear: gen ET2 = sum(servicepop*pop) / sum(pop)
by monthyear: replace ET2=ET2[_N]
//Average service population per clinic
drop if monthyear == 0
twoway line ET2 monthyear if tag, sort ///
xlabel(1 "1/2009" 13 "1/2010" 25 "1/2011" 37 "1/2012" 49 "1/2013" 61 "1/2014" 73 "1/2015" 85 "1/2016" 97 "1/2017"  ) ///
scheme(s1mono) name("graph_meanservicepop", replace) xtitle("Month-Year") ytitle("Mean Service Population Per Clinic (1,000s)") ///
xline(56 82, lcolor(red))title("Average Service Population")
graph export "Figure2_PanelB.png", replace


// FIGURE 3
// Map of change in service population

shp2dta using "${rawdata}/tl_2010_55_county10.shp", data("WI_data2")  coor("WI_coordinates2") genid(_ID) replace
use "${dtasave}/ClinicDataSet.dta", clear
 cap drop _merge
egen servicepop2009 = mean(servicepop_measure1/1000)if year == 2009, by(county) 
egen servicepop2016  = mean(servicepop_measure1/1000)if year == 2016, by(county)

preserve 
collapse (min) servicepop2016 servicepop2009, by (county)

gen changeservicepop = (servicepop2016- servicepop2009)/1000
tempfile temp1
 save "`temp1'" 
 restore
 merge m:1 county using "`temp1'" 
 preserve
 use "${rawdata}/WI_data", clear
 destring, replace
  g county = STATEFP10*1000+ COUNTYFP10
  
  tempfile temp2
   save "`temp2'" 
 restore
 drop _merge
  merge m:1 county using "`temp2'" 
  
preserve 
  import excel using "${rawdata}/WIregion2009_clinic_pointdata.xlsx", clear
  list


  rename A xcoord
  rename B ycoord
  rename C type

    save  temp4.dta, replace
  restore
  replace servicepop2009 = 179288  if id == 70 // Forest county is miscategorized
  preserve
  keep if year == 2009
  keep if month == 7
 spmap servicepop2009 using "${rawdata}/WI_coordinates", id(id) fcolor(Blues)   clmethod(custom) clbreaks(0 100 125 150 175 200 225 250 350 500)  ///
 legtitle("Women 15 to 44, per clinic") title("Average Service Population, 2009 ") name("WI_changeservicepopemap_2009", replace) ///
  point(data("temp4.dta") xcoord(ycoord)               ///
    ycoord(xcoord)     by(type)    ///
    psize(absolute) fcolor(green green*0.25 green)  ) ///
	polygon(data("${rawdata}/region2009_WI.dta") ocolor(black) osize(*2))
 *graph export "WI_changeservicepopemap_2009.png", replace
restore
  preserve
  keep if year == 2016
  keep if month == 7
 spmap servicepop2016 using "${rawdata}/WI_coordinates", id(id) fcolor(Blues)   clmethod(custom) clbreaks(0 100 125 150 175 200 225 250 350 500)  ///
 legtitle("Women 15 to 44, per clinic") title("Average Service Population, 2016 ") name("map2016", replace) ///
  point(data("temp4.dta") xcoord(ycoord)               ///
    ycoord(xcoord)    by(type)        ///
    psize(absolute) fcolor(green green red) ) ///
	polygon(data("${rawdata}/region2016_WI.dta") ocolor(black) osize(*2))
	
 *graph export "WI_changeservicepopemap_2016.png", replace
restore
 

 
// FIGURE 4
/*
comment in if using restricted data
do "WIAbortion_syntheticanalysis_birthrate"
*/
// FIGURE 5

do "${home}/WIabortion_syntheticanalysis.do"

// FIGURE 6
use "${home}/dtafiles/WI_distclinic_annual.dta", clear


sort county year
g distchange = closest_dist_[_n]-closest_dist_[_n-1] 
 g treatment_distincrease = (closest_dist_[_n]-closest_dist_[_n-1] > 25)

 tab year treatment_distincrease
 
xtset county year
 sort county year

cap drop treat_lag* treat_lead*

foreach num of numlist 1/8{
  g treat_lead`num'= treatment_distincrease[_n+`num']
   g treat_lag`num'= treatment_distincrease[_n-`num']
   }
 g treat_lag4plus = (treat_lag4 == 1| treat_lag5 == 1| treat_lag6 == 1| treat_lag7 == 1| treat_lag8 == 1)
g treat_lead4plus = (treat_lead4 == 1| treat_lead5 == 1| treat_lead6 == 1| treat_lead7 == 1| treat_lead8 == 1)

replace treat_lead1 = 0
xtset county year

xtpoisson abortion_county treat_lead4plus treat_lead3 treat_lead2 treat_lead1   treatment_distincrease treat_lag1 treat_lag2 treat_lag3 treat_lag4plus PCI unemp  pop*   i.year, fe i(county) vce(robust) exposure(pop)
estimates store leadslags_abortion

coefplot (leadslags_abortion) , keep(treat*) drop (treat_lag4plus) ///
order(treat_lead4plus treat_lead3 treat_lead2 . treatment_distincrease treat_lag1 treat_lag2 treat_lag3 ) ///
 xline(4.5, lpattern(dash) lcolor(gray))yline(0)lcolor(blue) ///
 graphregion(color(white)) bgcolor(white) ///
vertical label ///
 ytitle("Change in Abortion Rate") omitted  ///
nokey xtitle("Years Since Distance Change") xlabel(1 "t>=-4" 2 "t=-3" 3 "t=-2" 4 "t=-1" 5 "t=0" 6 "t=+1" 7 "t=+2" 8 "t=+3") title("Event Study: Increase in Distance > 25miles")  name("mile25_eventstudy", replace) 
graph save "Figure6.gph", replace
graph export "Figure6.png", replace

// FIGURE 7
/*
Comment in if you have access to NVSS

use WI_distclinic_notoffset.dta , clear
sort county monthyear
g distchange = closest_dist_[_n]-closest_dist_[_n-1] 
 g treatment_distincrease = (closest_dist_[_n]-closest_dist_[_n-1] > 25)

 tab monthyear treatment_distincrease
replace treatment_distincrease = . if monthyear <= 1
 drop if county == .
xtset county monthyear
 sort county monthyear

cap drop treat_lag* treat_lead*

foreach num of numlist 1/36{
  g treat_lead`num'= treatment_distincrease[_n+`num']
   g treat_lag`num'= treatment_distincrease[_n-`num']
   }
 g treat_lag25plus = (treat_lag25 == 1 | treat_lag26 == 1| treat_lag27 == 1| treat_lag28 == 1|  treat_lag29 == 1| treat_lag30 == 1| treat_lag31 == 1| treat_lag32 == 1|treat_lag33 == 1| treat_lag34 == 1| treat_lag35 == 1|treat_lag36 == 1 )
g treat_lead25plus = (treat_lead25 == 1 |treat_lead26 == 1| treat_lead27 == 1| treat_lead28 == 1|  treat_lead29 == 1| treat_lead30 == 1| treat_lead31 == 1| treat_lead32 == 1|treat_lead33 == 1| treat_lead34 == 1| treat_lead35 == 1|treat_lead36 == 1 )
xtset county monthyear
replace treat_lead1 = 0
xtset county monthyear

// Creating one-year moving average of birth rate and out of state birth %
sort county monthyear
foreach num of numlist 1/6{
bysort county: gen birthrate_lag`num' =birthrate[_n-`num']
bysort county: gen birthrate_lead`num' =birthrate[_n+`num']
}



xtpoisson allbirth treat_lead25plus treat_lead24 treat_lead23 treat_lead22 treat_lead21 treat_lead20 treat_lead19 treat_lead18 treat_lead17 treat_lead16 treat_lead15 treat_lead14 treat_lead13 treat_lead12 treat_lead11 treat_lead10 treat_lead9 treat_lead8 treat_lead7 treat_lead6 treat_lead5 treat_lead4 treat_lead3 treat_lead2 treat_lead1 treatment_distincrease treat_lag1 treat_lag2 treat_lag3 treat_lag4 treat_lag5 treat_lag6 treat_lag7 treat_lag8 treat_lag9 treat_lag10 treat_lag11 treat_lag12  treat_lag13  treat_lag14  treat_lag15  treat_lag16  treat_lag17  treat_lag18  treat_lag19  PCI unemp  pop* PPdist i.monthyear, fe i(county) vce(robust) exposure(pop)
estimates store leadslags_birth

coefplot (leadslags_birth, level(90)), keep(treat*) drop (treat*plus treat_lag13  treat_lag14  treat_lag15  treat_lag16  treat_lag17  treat_lag18  treat_lag19  treat_lag20  treat_lag21  treat_lag22  treat_lag23  treat_lag24) ///
order(treat_lead24 treat_lead23 treat_lead22 treat_lead21 treat_lead20 treat_lead19 treat_lead18 treat_lead17 treat_lead16 treat_lead15 treat_lead14 treat_lead13 treat_lead13 treat_lead12 treat_lead11 treat_lead10 treat_lead9 treat_lead8 treat_lead7 treat_lead6 treat_lead5 treat_lead4 treat_lead3 treat_lead2 . treatment_distincrease treat_lag1 treat_lag2 treat_lag3 treat_lag4 treat_lag5 treat_lag6 treat_lag7 treat_lag8 treat_lag9 treat_lag10 treat_lag11 treat_lag12  ) ///
 yline(0)  lcolor(blue) ///
 graphregion(color(white)) bgcolor(white) ///
vertical label ///
xlabel(1 "t= -24" 7 "t= -18" 13 "t= -12" 19 "t=-6" 25 "t= 0" 31 "t=+6" 37 "t=+12" ) xtitle("Months Since Event") xline(24.5, lpattern(dash) lcolor(gray)) ///
 ytitle("Change in  Birth Rate") omitted  ///
nokey  title("Event Study: Increase in Distance >25 miles")  name("mile25_eventstudy_birth", replace) 
graph save "Figure7.gph", replace
graph export "Figure7.png", replace

*/
// FIGURE 8
use "${home}/dtafiles/WI_distclinic_annual.dta", clear


capture program drop _all
program marginsp
args access xll xd xup delta name
# delimit ;
margins, 
expression(
	100*(
	exp(
		_b[`access']*((`access'+`delta')-`access')+
		_b[`access'SQ]*((`access'+`delta')^2-`access'^2)
		)
	-1)
)
at(`access'=(`xll'(`xd')`xup'))
saving(`name', replace)
 ;
# delimit cr
end


g travel_distance = closest_dist_ 
gen travel_distanceSQ=travel_distance*travel_distance
gen travel_distanceCU=travel_distance*travel_distance*travel_distance
gen travel_distanceQU=travel_distance*travel_distance*travel_distance*travel_distance
xtset countyfp year
eststo: xtpoisson abortion_county travel_distance travel_distanceSQ PCI unemp  PPdist pop* i.year,  fe  i(countyfp) vce(robust) exposure(pop)
local d = 50
marginsp travel_distance 0 25 150 `d' td
marginsplot,  ytitle("% Effect of a `d'-Mile Increase in Distance") xlabel(0 "0" 50 "50" 100 "100" 150 "150" )  xtitle("Base Travel distance to nearest clinic") title("") yline(0) legend(pos(6)) name("abortion", replace) ///
plotregion(fcolor(white)) graphregion(fcolor(white))  title("Change in Abortion Rate, by Starting Distance")
graph export "${WIdata}/Table8.png", replace
	
 
// FIGURE 9 -- Uncomment if using restricted data
/*
use "WI_distclinic.dta", clear
   
capture program drop _all
program marginsp
args access xll xd xup delta name
# delimit ;
margins, 
expression(
	(
	exp(
		_b[`access']*((`access'+`delta')-`access')+
		_b[`access'SQ]*((`access'+`delta')^2-`access'^2)
						)
	-1)
)
at(`access'=(`xll'(`xd')`xup'))
saving(`name', replace)
 ;
# delimit cr
end

g travel_distance = closest_dist_
gen travel_distanceSQ=travel_distance*travel_distance
gen travel_distanceCU=travel_distance*travel_distance*travel_distance
gen travel_distanceQU=travel_distance*travel_distance*travel_distance*travel_distance
gen travel_distanceQUI=travel_distance*travel_distance*travel_distance*travel_distance*travel_distance
xtset county monthyear
eststo: xtpoisson allbirth travel_distance travel_distanceSQ  PCI unemp pop* PPopen,  fe  i(county) vce(robust) exposure(pop)
local d = 50
marginsp travel_distance 0 25 125 `d' td
marginsplot,  ytitle("% Effect of a `d'-Mile Increase in Distance") xlabel(0 "0" 50 "50" 100 "100" 150 "150" )  xtitle("Travel distance to nearest clinic") title("") yline(0) legend(pos(6)) name("Birth", replace) ///
plotregion(fcolor(white)) graphregion(fcolor(white)) title("Change in Birth Rate, by Starting Distance")
graph export "${WIdata}/Figure9.png", replace
*/
// FIGURE 10
import excel using "${rawdata}/OOSabortiondata.xlsx", firstrow clear
label var WIRes_IL "Abortions to WI res. in IL"
label var WIRes_MN "Abortions to WI res. in MN"
destring, replace force
sort Year
twoway (line WIRes_IL Year, lcolor(orange) lwidth(medthick)) (line WIRes_MN Year, lcolor("0 158 115") lpattern(dash) lwidth(medthick)), ///
 name("graph_outofstate", replace) xtitle("Year") ytitle("Count of Abortions") ///
xline(2013 2015.5, lcolor(grey)) title("Number of Abortions to WI residents in Neighboring States") ///
plotregion(fcolor(white)) graphregion(fcolor(white)) 
graph export "Figure10.png", replace


// TABLE 1
local restrict = 0 // set equal to 1 if using data with births; set equal to 0 if only using public use data
if `restrict' == 0 {
use "${home}/dtafiles/WI_distclinic_annual.dta" , clear
}
if `restrict' == 1{
	use "V:\FletcherBirthRecords\Joanna\Wisconsin\WI_distclinic_annual.dta" , clear
}
cap drop Perc_White
 cap drop Perc_Black
 cap drop Perc_OtherRace
 g Perc_White = popWhiteWomen/(pop)
 g Perc_Black = popBlackWomen/(pop)
 g Perc_OtherRace = 1- Perc_White - Perc_Black
 replace pop = pop_F15to19 + pop_F20to24 + pop_F25to29 + pop_F30to34+ pop_F35to39 + pop_F40to44
 cap drop Perc_1* Perc_2* 
 cap drop Perc_3*
 g Perc_15to19 = pop_F15to19/pop
 g Perc_20to24 = pop_F20to24/pop
 g Perc_25to29 = pop_F25to29/pop
 g Perc_30to34 = pop_F30to34/pop
 g Perc_35to39 = pop_F35to39/pop
 g Perc_40to44 = pop_F40to44/pop
 if `restrict' ==1 {
 // Birthrate by race and age group
 cap drop BR_*
g BR_white= 1000*whitebirth/popWhiteWomen
g BR_black= 1000*blackbirth/popBlackWomen
g BR_teen = 1000*teenbirth/pop_F15to19
g BR_twenties = 1000*twentybirth/(pop_F20to24+pop_F25to29)
g BR_thirties = 1000*thirtyplusbirth/(pop_F30to34+pop_F35to39 +pop_F40to44)
 }
 g  servicepop_measure=  servicepop_measure1/1000
global summstat_var AbortionRate closest_dist_ servicepop_measure ///
Perc_White Perc_Black Perc_OtherRace Perc_15to19 Perc_20to24 Perc_25to29 Perc_30to34 ///
Perc_35to39 Perc_40to44  ///
unemp PCI
global summstat_restrict birthrate BR_white BR_black BR_teen BR_twenties BR_thirties
eststo clear
svyset county [pw=pop]
// Weighted by number of women in each county
if `restrict' == 0{
estpost sum $summstat_var if year < 2017 [w=pop], detail
eststo sum_whole
estpost sum $summstat_var if year == 2012 [w=pop], detail
eststo sum_2012
estpost sum $summstat_var if year == 2016 [w=pop], detail
eststo sum_2016
}
if `restrict' == 1{
estpost sum $summstat_var $summstat_restrict if year < 2017 [w=pop], detail
eststo sum_whole
estpost sum $summstat_var $summstat_restrict if year == 2012 [w=pop], detail
eststo sum_2012
estpost sum $summstat_var $summstat_restrict if year == 2016 [w=pop], detail
eststo sum_2016
}

label var AbortionRate "Abortion Rate, per 1000 women"
label var closest_dist_ "Distance to Closest Clinic"
label var servicepop_measure "Average Service Pop. (1000s)"
label var Perc_White "Women 15 to 44, Percent White"
label var Perc_Black "Women 15 to 44, Percent Black"
label var Perc_OtherRace "Women 15 to 44, Percent Other Race"
label var Perc_15to19 "Percent Age 15 to 19"
label var Perc_20to24 "Percent Age 20 to 24"
label var Perc_25to29 "Percent Age 25 to 29"
label var Perc_30to34 "Percent Age 30 to 34"
label var Perc_35to39 "Percent Age 35 to 39"
label var Perc_40to44 "Percent Age 40 to 44"
if `restrict' == 1{
label var birthrate "Birth Rate, per 1000 women"
label var BR_white "Birth Rate, White"
label var BR_black "Birth Rate, Black"
label var BR_teen "Birth Rate, 15 to 19"
label var BR_twenties "Birth Rate, 20 to 29"
label var BR_thirties "Birth Rate, 30 and over"
}
label var unemp "County Unemployment Rate"
label var PCI "County Per Capita Income"


		esttab sum_whole sum_2012 sum_2016  ///
		 using "${WIdata}/Table1.tex", ///
		cell("mean(fmt(3)) sd(fmt(3) par) ") nonumbers ///
		mtitles("2009-2016" "2012" "2016" ) unstack ///
		label title("Summary Statistics") replace
		
		
// TABLE 2
use "${home}/dtafiles/WI_distclinic_annual.dta" , clear
 eststo clear
 xtset, clear
 xtset county year
 cap drop dist dist_sq
 g dist= closest_dist_/100
 g dist_sq = dist^2
eststo: xtpoisson abortion_county dist PCI unemp pop* PPdist i.year if abortion_county != 3,  fe  i(county) vce(robust) exposure(pop)
eststo: xtpoisson abortion_county dist dist_sq PCI unemp pop* PPdist i.year if abortion_county != 3,  fe  i(county) vce(robust) exposure(pop)
 eststo: xtpoisson abortion_county Clinic_25to50 Clinic_50to100 NoClinic_0to100 PCI unemp PPdist pop* i.year if abortion_county != 3, fe i(county) vce(robust) exposure(pop)
  esttab using "${WIdata}/Table2.tex",  replace se 
// TABLE 3 Uncomment if you have access to restricted data
/*

use "WI_distclinic.dta", clear
eststo clear
drop if county == .
xtset county monthyear
 cap drop dist dist_sq
 g dist= closest_dist_/100
 g dist_sq = dist^2
eststo: xtpoisson allbirth dist PCI unemp pop_F15to19 pop_F20to24 pop_F25to29 popWhite popBlack popNat* PPdist i.year ,  fe  i(county) vce(robust) exposure(pop)
eststo: xtpoisson allbirth dist dist_sq PCI unemp pop_F15to19 pop_F20to24 pop_F25to29 popWhite popBlack popNat*  PPdist i.year,  fe  i(county) vce(robust) exposure(pop)
eststo: xtpoisson allbirth Clinic_25to50 Clinic_50to100 NoClinic_0to100 PCI unemp pop_F15to19 pop_F20to24 pop_F25to29 popWhite popBlack popNat*  PPdist i.year, fe  i(county) vce(robust) exposure(pop)

 esttab using Table3.tex,  replace se 
 eststo clear
*/

// TABLE 4
** Col 1 and 2:
use "${home}/dtafiles/WI_distclinic_annual.dta" , clear
 cap drop servicepop_1000
g servicepop_1000 = servicepop_measure1/1000 
 eststo clear
xtset county year
 cap drop dist dist_sq
 g dist= closest_dist_/100
 g dist_sq = dist^2
eststo: xtreg AbortionRate  servicepop_1000  PCI unemp pop* PPdist i.year, fe   vce(cluster county)
eststo: xtreg AbortionRate dist   servicepop_1000  PCI unemp pop* PPdist i.year, fe   vce(cluster county)
eststo: xtpoisson AbortionRate  servicepop_1000  PCI unemp pop* PPdist i.year,  fe  i(county) vce(robust) exposure(pop)

eststo: xtpoisson AbortionRate dist  servicepop_1000  PCI unemp pop*  PPdist i.year, fe i(county) vce(robust) exposure(pop)
 esttab using "${WIdata}/Table4_Col1and2.tex",  replace se 

 
// TABLE 5
/* Comment in if you have access to restricted data
use "WI_distclinic.dta", clear
 g BR_teen = teenbirth/pop_F15to19
 g BR_twenties = twentybirth/(pop_F20to24 + pop_F25to29)
 g BR_thirtyplus = thirtyplusbirth/(pop_F30to34 + pop_F35to39 + pop_F40to44)
 g BR_black = blackbirth/ popBlack
 g BR_white = whitebirth/ popWhite
 g nonwhitebirth = blackbirth +hispbirth+ asianbirth
g dist = closest_dist_/100
 eststo clear
xtset county monthyear
foreach var of varlist teenbirth twentybirth thirtyplusbirth whitebirth nonwhitebirth unmarriedbirth marriedbirth{
eststo: xtpoisson `var' dist PCI unemp pop* PPdist i.monthyear,  fe  i(county) vce(robust) exposure(pop)

}
 
 esttab using Table5.tex,  replace se 
*/
// TABLE 6
use "${home}/dtafiles/WI_distclinic_annual.dta" , clear
g dist = closest_dist_/100
eststo clear
 eststo: xtpoisson abortion_county  c.closest_dist_##Act217Post PCI unemp pop* i.year if abortion_county != 3 ,  fe  i(county) vce(robust) exposure(pop)
  eststo: xtpoisson abortion_county Act217Post##(Clinic_25to50 Clinic_50to100 NoClinic_0to100) PCI unemp pop* i.year if abortion_county != 3, noomit fe i(county) vce(robust) exposure(pop)
  esttab using "${WIdata}/Table6.tex",  replace se star(+ 0.1 * 0.05 ** 0.01 ) 

**** APPENDIX ****

//Figure A1
use "${home}/dtafiles/WI_distclinic_annual.dta" , clear
 eststo clear
 xtset, clear
 xtset county year
capture program drop _all
program marginsp
args access xll xd xup delta name
# delimit ;
margins, 
expression(
	100*(
	exp(
		_b[`access']*((`access'+`delta')-`access')+
		_b[`access'SQ]*((`access'+`delta')^2-`access'^2)
		)
	-1)
)
at(`access'=(`xll'(`xd')`xup'))
saving(`name', replace)
 ;
# delimit cr
end

 eststo: xtpoisson abortion_county Clinic_25to50 Clinic_50to100 NoClinic_0to100 PCI unemp  PPdist pop* i.year  if abortion_county != 3, fe i(county) vce(robust) exposure(pop)
 coefplot, keep(Clinic_25to50 Clinic_50to100 NoClinic_0to100) graphregion(color(white)) bgcolor(white) vertical yline(0) yscale(range (-0.5 0.01)) label ///
 ytitle("Percent Change in Annual Abortions") title("Change in County Abortion Rates by Distance Bin") name(distancebin, replace)
 graph export "${WIdata}/FigureA1.png", replace
// Figure A2

//Figure A3
import excel using ${rawdata}/OOSabortiondata.xlsx, firstrow clear
label var MIRes_WI "Abortions to MI res. in WI"
label var UP_MI  "Abortions to UP res. in MI"
destring, replace force
sort Year
twoway (line MIRes_WI Year, lcolor(orange) lwidth(medthick)) (line UP_MI Year, lcolor("0 158 115") lpattern(dash) lwidth(medthick) ) if Year < 2018, ///
 name("graph_Michigan", replace) xtitle("Year") ytitle("Count of Abortions") ///
xline(2013 2015.5, lcolor(grey)) title("Number of Abortions to Michigan: In WI and In MI") ///
plotregion(fcolor(white)) graphregion(fcolor(white)) 
graph export "FigureA3.png", replace

// Table A-1
/*
Only runs if you have restricted access NVSS files
do "${home}/StateComparison_Demographics.do"

*/

// Table A-2
use "${home}/dtafiles/WI_distclinic_annual.dta" , clear
cap drop dist_* abort_*
sort county year
egen dist_start = mean(closest_dist_)if year == 2009, by(county) 
egen dist_pre2013close  = mean(closest_dist_)if year == 2013, by(county)
egen dist_end = mean(closest_dist_)if year == 2017, by(county)
egen abort_start = mean(AbortionRate)if year == 2009, by(county) 
egen abort_pre2013close  = mean(AbortionRate)if year ==2013, by(county)
egen abort_end = mean(AbortionRate)if year == 2017, by(county)
eststo clear
preserve 
collapse (min) dist_start dist_pre2013close dist_end abort_start abort_pre2013close abort_end, by (county)

gen dist0 = dist_pre2013close-dist_start
gen dist1 = dist_end - dist_pre2013close
gen abort0 = abort_pre2013close-abort_start
gen abort1 = abort_end - abort_pre2013close
eststo: reg dist1 abort0
 esttab using "${WIdata}/TableA3_col1.tex",  replace se 


 restore
 * Columm 2 requires restricted access data; uncomment if have access
 /*
 use "${home}/dtafiles/WI_distclinic.dta" , clear
cap drop dist_* birth_*
sort county monthyear
egen dist_start = mean(closest_dist_)if monthyear == 1, by(county) 
egen dist_pre2013close  = mean(closest_dist_)if monthyear == 55, by(county)
egen dist_end = mean(closest_dist_)if monthyear == 102, by(county)
egen birth_start = mean(birthrate)if monthyear == 1, by(county) 
egen birth_pre2013close  = mean(birthrate)if monthyear == 55, by(county)
egen birth_end = mean(birthrate)if monthyear == 102, by(county)

preserve 
collapse (min) dist_start dist_pre2013close dist_end birth_start birth_pre2013close birth_end, by (county)

gen dist0 = dist_pre2013close-dist_start
gen dist1 = dist_end - dist_pre2013close
gen birth0 = birth_pre2013close-birth_start
gen birth1 = birth_end - birth_pre2013close
eststo: reg dist1 birth0

 restore
 
 esttab using TableA2_col2.tex,  replace se 
 eststo clear
 */
// Table A-3
use "${home}/dtafiles/WI_distclinic_annual.dta" , clear
 eststo clear
 xtset, clear
xtset county year
eststo: xtreg AbortionRate closest_dist_ PCI unemp pop* PPdist i.year if abortion_county != 3,  fe   vce(cluster county)
eststo: xtreg AbortionRate closest_dist_ distsq PCI unemp pop* PPdist i.year if abortion_county != 3,  fe   vce(cluster county)
eststo:xtreg AbortionRate Clinic_25to50 Clinic_50to100 NoClinic_0to100 PCI unemp pop* PPdist i.year if abortion_county != 3 , fe   vce(cluster county)

 esttab using "${WIdata}/TableA3.tex",  replace se 
// Table A-4
/* Uncomment if have access to NVSS data

use "WI_distclinic.dta", clear
eststo clear
drop if county == .
xtset county monthyear
 cap drop dist dist_sq
 g dist= closest_dist_/100
 g dist_sq = dist^2
eststo: xtreg allbirth dist PCI unemp pop_F15to19 pop_F20to24 pop_F25to29 popWhite popBlack popNat* PPdist i.year ,  fe   vce(cluster county)
eststo: xtreg allbirth dist dist_sq PCI unemp pop_F15to19 pop_F20to24 pop_F25to29 popWhite popBlack popNat*  PPdist i.year,  fe   vce(cluster county)
eststo: xtreg  allbirth Clinic_25to50 Clinic_50to100 NoClinic_0to100 PCI unemp pop_F15to19 pop_F20to24 pop_F25to29 popWhite popBlack popNat*  PPdist i.year,  fe   vce(cluster county)

 esttab using Table3.tex,  replace se 
 eststo clear
*/
// Table A-5
/* Uncomment if have access to NVSS data
use "WI_distclinic_9monthsoffset.dta", clear

eststo clear

drop if county == .
xtset county monthyear
eststo: xtpoisson allbirth closest_dist_ PCI unemp pop_F15to19 pop_F20to24 pop_F25to29 popWhite popBlack popNat* PPdist i.monthyear ,  fe  i(county) vce(robust) exposure(pop)
use "WI_distclinic_12monthsoffset.dta", clear

drop if county == .
xtset county monthyear
eststo: xtpoisson allbirth closest_dist_ PCI unemp pop_F15to19 pop_F20to24 pop_F25to29 popWhite popBlack popNat* PPdist i.monthyear ,  fe  i(county) vce(robust) exposure(pop)

esttab TableA5.tex,  replace se ///
star(+ 0.10 * 0.05 ** 0.01) keep(closest_dist_)
*/
//Table A-6
use "${home}/dtafiles/WI_distclinic_annual.dta" , clear

egen dist2009 = mean(closest_dist_)if year == 2009, by(county) 
egen dist2017  = mean(closest_dist_)if year == 2017, by(county)
drop _merge 
cap drop changedist
preserve 
collapse (min) dist2017 dist2009, by (county)

gen changedist = dist2017- dist2009
tempfile temp1
 save "`temp1'" 
 restore
 merge m:1 county using "`temp1'" 

g distchange_group = 1 if changedist <=0

xtile distchangeterc = changedist if changedist > 0, nq(3)
replace distchange_group = distchangeterc+1 if changedist > 0

 xtset county year
eststo: xtpoisson abortion_county closest_dist_ PCI unemp pop* PPdist i.year if abortion_county != 3 & distchange_group != 3,  fe  i(county) vce(robust) exposure(pop)

 eststo: xtpoisson abortion_county Clinic_25to50 Clinic_50to100 NoClinic_0to100 PCI unemp pop* PPdist i.year  if abortion_county != 3 & distchange_group !=3,  fe i(county) vce(robust) exposure(pop)
 
 use "${home}/dtafiles/WI_distclinic_annual.dta" , clear
gen NotOutOfState= (ClinicID < 7)
bysort county: egen NeverOOS = min(NotOutOfState)


 xtset county year
eststo: xtpoisson abortion_county closest_dist_ PCI unemp pop* PPdist i.year if abortion_county != 3 & NeverOOS ==1,  fe  i(county) vce(robust) exposure(pop)
eststo: xtpoisson abortion_county closest_dist_ distsq PCI unemp pop* PPdist i.year  if abortion_county != 3& NeverOOS ==1,   fe  i(county) vce(robust) exposure(pop)
 eststo: xtpoisson abortion_county Clinic_25to50 Clinic_50to100 NoClinic_0to100 PCI unemp pop* PPdist i.year  if abortion_county != 3& NeverOOS ==1,  fe i(county) vce(robust) exposure(pop)

*esttab using TableA6.tex,  replace se ///
*star(+ 0.10 * 0.05 ** 0.01)  keep(closest_dist_)
